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SOLUTION  OF  TRANSCENDENTAL  AND  ALGEBRAIC  EQUATIONS 
WITH  APPLICATION  TO  WAVE  PROPAGATION  IN  ELASTIC  PLATES 

INTRODUCTION 

Transcendental  equations  with  the  general  form 

A(x)  =  C(x)  (1) 

B(x)  D(x) 

arise  in  elastic  plate  theory  as  the  characteristic  equation  governing  the 
propagation  of  symmetric  and  antisymmetric  waves.  This  report  documents 
the  development  of  a  generalized  iterative  rootfinder  for  the  solution  of 
this  type  of  equation  and  its  application  to  elastic  plate  theory.  This 
report,  therefore,  represents  a  dichotomy  of  topics.  The  initial  portion 
of  the  text  gives  a  thorough  development  of  the  generalized  rootfinder  used  in 
solving  transcendental  equations  in  the  form  of  Eq.  (1).  The  latter  part 
applies  the  generalized  method  to  solve  a  number  of  equations  arising  in 
elastic  plate  theory. 

The  organization  of  this  paper  is  as  follows:  First,  the  approximation 
of  a  root  of  the  equation  f(x)  ■  0  by  a  known  iterative  method  is  reviewed. 

A  simple  transcendental  equation  is  then  chosen  to  exemplify  this  basic 
method.  Next,  the  generalized  iterative  rootfinder  is  introduced.  From  ex¬ 
amining  this  generalized  method,  the  reader  will  see  how  the  basic  iterative 
scheme  reviewed  previously  can  be  successively  applied  to  determine  the  roots 
of  the  general  equation  with  the  form  of  Eq.  (1).  Using  this  general  itera¬ 
tive  rootfinder,  one  can  calculate  and  plot  dispersion  curves  for  both  the 
Manuscript  submitted  May  22,  1981. 


flexural  and  extensional  waves  In  rubber,  steel,  and  beryllium  plates.  These 
curves  are  found  to  approach  a  real  value  <  (ratio  of  the  Rayleigh  to  shea” 
wave  speed)  for  the  two  fundamental  inodes  of  plate  vibration.  The  values 
of  k  calculated  from  exact  elasticity  theory  are  compared  with  those 
obtained  from  an  approximate  formula  by  Victorov.  Tn  addition,  the  inflection 
and  maximum  points  of  the  strain  distribution  throughout  a  free  plate  were 
calculated.  Finally,  the  strain  distribution  throughout  a  free  plate  was 
determined  and  results  presented. 

ITERATIVE  ROOTFINDER 

Basic  Method 

Suppose  F(x)  is  a  real  continuous  function  such  that 

F(x)  =  0  (2) 

has  real  roots.  In  general,  these  real  roots  can  be  computed  to  any  degree 
of  accuracy  by  an  iterative  procedure.  The  Iterative  method  used  is  based 
on  rewriting  Eq.  (2)  in  the  form 

f(x)  =  g(x).  (3) 

Let  z  be  an  actual  root  of  Eq.  (3).  Take  x^  as  the  initial  approximation 
of  z.  (This  initial  trial  root  x^  will  subsequently  be  referred  to  as  the 
"seed"  for  the  iterative  calculation.)  If  x  is  replaced  by  x^  In  the  right- 
hand  side  of  Eq.  (3),  then 

f(x)  =  g(xQ),  (4) 

which  can  be  solved  by  calculating  x^  in  th>'  form 

Xj  -  f  1Cg(x0)].  (5) 
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The  value  x^  obtained  by  Eq.  (5)  becomes  the  argument  of  a  new  approximation 

x2  =  f  1  [ g(x1 ) ] ,  (6) 

♦*  Vi 

where  the  value  x  that  satisfies  Eq.  (6)  is  x^>  In  general,  the  nu 
approximation  x^  is  determined  by  solving 

x  =  f  1[g(x  )].  (7a) 

n  n-i 

If  the  sequence  x^,  x  ,  x^  •••  tends  to  a  limit  z  and  if  f(x)  and  g(x)  are 
continuous  at  the  point  z,  then  z  is  the  root  of  the  given  equation. 

One  might  suppose  that  Eq.  (5)  could  have  been  written  as 

xL  =  g  1[f(xQ)],  (7b) 

instead  of  Eq.  (7a).  However,  the  side  of  Eq.  (3)  that  needs  to  remain 
fixed  by  the  insertion  of  x^  is  not  arbitrary.  If  the  two  functions  f(x) 
and  g(x)  have  first  derivatives  such  that 

| f * (x) |  >  |g’(x) | ,  (8) 

then  the  sequence  Xq,  x^,  x^  ...  is  convergent  whenever  written  in  the  form 
of  Eq.  (7a),  but  not  when  written  in  the  form  of  Eq.  (7b).  That  is,  the 
function  with  the  smaller  slope  in  the  neighborhood  of  z  of  Eq.  (3)  must 
be  fixed  in  order  to  obtain  successively  better  approximations  of  the  true 
root  z. 

Sokolnikoff  and  Redheffer  [1]  describe  the  geometric  interpretation  of 
this  iterative  procedure.  Equation  (3)  describes  two  curves,  namely 
y  ■  f(x)  and  y  ■  g(x),  which  intersect  at  the  root  z.  If  the  slopes  of 
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these  two  curves  have  the  same  sign  and  satisfy  Eq.  (8),  the  convergence  is 
in  a  step-like  pattern  as  shown  in  Fig.  1.  One  may  note  that  in  this  case 


Fig.  1.  Intersection  of  two  equations  having  slopes  with  the  same  sign. 

Solid  line  is  result  of  fixing  equation  with  smaller  slope. 

Dashed  line  indicates  a  convergence  as  a  result  of  choosing 
the  larger  slope. 


successive  approximations  approach  the  true  root  z  from  one  side.  On  the 
other  hand,  when  the  slopes  of  the  curves  have  opposite  signs,  as  shown  in 


Fig.  2,  convergence  has  a  spiral-type  configuration.  The  successive  approx¬ 
imations  in  Fig.  2  oscillate  on  both  sides  of  z  and  yet  approach  the  true 
root  z  when  the  curve  g(x),  with  the  smaller  slope  is  fixed.  The  dotted  line 
shows  how  the  iterative  process  can  generate  approximations  that  diverge  from 
z.  This  occurs  when  f(x),  the  curve  with  the  larger  slope,  is  fixed  rather 
than  g(x).  Regardless  of  whether  convergence  takes  place  in  a  stepwise  or  a 


Fig.  2.  Intersection  of  two  equations  having  slopes  of  opposite  sign. 

spiral  fashion,  the  curve  with  the  smaller  slope  must  be  the  one  fixed  in 
Eq.  (3),  otherwise  consecutive  approximations  will  not  approach  z  but  instead 
diverge  from  it. 

Since  with  any  numerical  method,  exact  roots  cannot  be  found,  a  degree 
of  accuracy  must  be  decided  on.  Therefore,  whenever  the  difference  lxn  “  xn_ 


Is  within  a  predetermined  tolerance,  is  accepted  as  a  sufficiently  good 
approximation  of  z,  the  true  root  of  Eq.  (2). 


Example 

It  seems  appropriate  to  examine  at  this  time  an  example  of  this  iter¬ 
ative  scheme  [2],  Suppose  one  wishes  to  determine  the  roots  of 

tan(x)  -  —  =  0.  (9) 

x 


The  real  roots  of  this  equation  are  the  absissas  of  the  point  of  inter¬ 
section  of  the  curves 


and 


y  =  tan (x) 


(10a) 

(10b) 


A  trial  root  or  seed  Xq  should  be  in  the  vicinity  of  the  true  root  for  a 
more  rapid  convergence.  However,  with  the  use  of  the  computer,  the  extra 
time  required  when  x^  is  not  chosen  near  the  actual  root  z  is  usually 
negligible. 


The  first  trial  root  Xq  may  be  chosen  to  be  equal  to  one.  The  slope 

2 

of  y  =  tan  x  is  greater  than  that  of  y  ■  -  at  x  a  1.0.  Therefore,  Eq.  (9) 
is  written  in  the  form 

x.^  *  tan  1(2/xq)  (11a) 

or  in  general,  for  the  n^  iteration 


x 

n 


( lib) 


The  first  seven  terms  of  the  sequence  generated  by  Eq.  (lib)  with  x^  *  1.0, 
are  x^  -  1.1071,  x2  *  1.0652,  x3  -  1.0814,  x^  -  1.0751,  x$  -  1.0775, 
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■» 


x.  =  1.0766,  and  x7  =  1.0769.  The  computations  can  he  terminated  at  this 
6  ' 

-k 

stage  if  the  predetermined  tolerance  t  is  set,  say,  at  3  *  10  because 

lx  -  x  |  <  t.  It  is  important  to  point  out  that  the  tolerance  is  not 
1  6  5 1  — 

t  Vl 

necessarily  the  error  between  the  n  computed  root  x^  and  the  true  root 
z,  as  indicated  in  many  books  (see  Appendix  A  for  a  detailed  analysis  of 
the  error  and  tolerance  size). 


It  should  also  be  noted  that  if  one  function  in  F,q.  (3)  has  a  larger 
slope  at  one  point,  this  does  not  ensure  that  it  will  have  the  larger  slope 

at  every  point.  For  instance,  in  the  example  considered  here,  the  deriva- 

2  -2 
tive  of  2/x  is  -2/x  and  the  derivative  of  tan(x)  is  [cos(x)J  .  When 

x  =  1,  the  slope  of  2/x  is  -2  while  the  slope  of  tan(x)  is  3.425.  Therefore, 


—  tan(x)  > 
dx 


dx  \x  J 


at  x  =  1.  However,  if  instead  0.1  were  chosen  as  the  value 


of  x,  the  slope  of  2/x  is  -200  while  the  slope  of  tan(x)  is  1.005.  Thus, 


one  would  find  that 


dx 


(!) 


>  —  tan(x)  at  x  =*  0.1.  Hence  by  Eq.  (8),  Eq.  (9) 
dx 


should  be  written  as 


x^  =  2  tanCxg)  (11c) 

rather  than  in  the  form  of  Eq.  (11c).  However,  Eq.  (11c)  will  produce 
successive  roots  that  diverge  from  the  true  root.  Therefore,  it  is  important 
that  the  initial  trial  root  Xq  be  close  enough  to  the  true  root  so  that  the 
slopes  of  the  function  f  and  g  in  Eq.  (3)  in  the  neighborhood  of  the  trial 
root  will  have  the  same  relation  to  one  another  as  they  have  in  the  neighbor¬ 


hood  of  the  true  root. 


Generalized  Method 


The  basic  method  discussed  above  is  effective  when  an  equation  can  be 
written  in  the  form  f(x)  =  g(x).  However,  for  an  equation  written  in  the 
form 

A(x )  _  C(x )  (12a) 

B(x)  D(x)  ’ 

the  iterative  rootfinder  described  above  must  be  made  more  general.  In  the 
generalized  method,  the  basic  iterative  method  described  above  is  successively 
applied  to  combinations  of  the  functions  that  appear  in  Eq.  (12)  until 
the  root  is  determined.  The  rationale  for  solving  Eq.  (12a)  is  as  follows: 

A  (x  )  C  (x  ) 

Let  =  F(x)  and  |^ry  =  G(x),  and  suppose  the  slopes  of  F(x)  and 

G(x)  can  be  determined.  If  Xq  is  the  first  trial  root  and  F(x)  has  the 
greater  slope  at  Xq,  then  Eq.  (12a)  can  be  written  in  the  form 

rey  ’  C<V-  <12h) 

Equation  (12b)  can  also  be  written  as 

A(x)  =  G(xq)B(x),  (13a) 

where  G(xq)  serves  as  a  constant  since  G  is  evaluated  at  x^.  Note  that  if 
G(x)  had  a  slope  larger  than  F(x)  at  Xq,  Eq.  (12a)  would  have  been  written 
as 

C(x)  =  F(xq)D(x).  (13b) 

The  two  pairs  of  functions,  y  =  A(x)  and  y  =  G(xq)B(x)  of  Eq.  (12a),  or 
y  =  C(x)  and  y  =  F(xq)D(x)  of  Eq.  (13b),  also  demand  a  slope  comparison  to 

assure  that  successive  approximations  will  approach  the  true  root  z.  Conse¬ 
quently  for  an  equation  in  the  form  of  (12a),  two  slope  comparisons  need 


to  be  made  for  each  approximation.  To  deal  with  an  equation  in  this  form, 
a  "Little  Loop-Big  Loop  (LL-BL)  Procedure"  was  designed.  The  "Little  Loop- 
Big  Loop  Procedure"  breaks  down  Eq.  (12a)  into  the  form  of  Eq.  (3). 
Basically  this  is  done  by  making  two  slope  comparisons  and  implementing  the 
iterative  procedure. 


It  should  be  kept  in  mind  that  this  LL-BL  procedure  works  not  only  for 

equations  of  the  form  of  Eq.  (12a),  but  can  also  be  implemented  for  equations 

A(x) 

involving  three  functions,  ^x)  =  C(x)»  tw0  functions  A(x)  =  B(x),  or  even 
one  function  such  as  A(x)  =  C  where  C  is  a  constant. 


First  the  slopes  of  F(x)  and  G(x)  are  compared  in  the  BL.  Then  by  the 
criterion  given  in  Eq.  (8),  Eq.  (12a)  will  be  written  in  the  form  of  either 
Eq.  (13a)  or  Eq.  (13b)  and  will  be  passed  to  the  LL.  In  the  LL,  the  slopes 
of  the  functions  of  Eq.  (13)  will  be  calculated.  Again  by  the  criterion 
given  in  Eq.  (8),  one  of  the  following  equations  will  result.  If  Eq.  (13a) 
was  used,  then  one  will  have 

x  =  A  1[B(x)  G(xq)  ]  (lAa) 

when 

|[G(x0)B(x)J'|  <  |A(x)'| 

or 

x  =  B  1[A(x)/G(xf)) )  (l4b) 

when 

| [G(xQ)B(x)  ]' j  >  5 A(x) '  | , 
and,  if  Eq.  (12b)  were  used, 

x  =  C  ^(F(xq)D(x)  ]  (14c) 
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when 

| (F<x0)D(x)]*  |  •  |c(x) '  | 
or 

x  =  D_1[C(x)/F(x0)]  (I4d) 

when 

|[F(x0)D(x)J*  |  >  | C (x ) '  |  . 

As  seen  above,  there  are  four  possible  ways  to  rewrite  F.q.  (12a);  however, 
there  is  only  one  equation  that  results  in  convergent  approximations  as 
discussed  earlier  in  the  basic  iterative  procedure. 

After  the  proper  version  of  Eq.  (14)  is  chosen,  the  difference  between 
successive  approximations  must  be  made  less  than  a  preset  tolerance.  This 
tolerance  is  called  the  LL-tolerance.  The  root  r  that  satisfies  the  LL 
tolerance  is  then  recirculated  into  Eq.  (12a)  to  again  make  the  BL  slope 
comparison.  After  the  proper  forms  of  Eqs.  (13)  and  (14)  have  been  chosen 
for  trial  root  r^,  the  iterative  procedure  is  carried  out  and  a  new  root  r2 
results,  where  r^  replaces  Xq.  In  general,  after  n  such  "recirculations", 
when  |r  .  -  r  I  is  less  than  a  preset  tolerance,  r  is  taken  to  be  the  root 
of  Eq.  (12a)  and  this  tolerance  is  called  the  BL  tolerance.  Note  there  is 
both  a  LL  tolerance  and  a  BL  tolerance,  and  these  need  not  be  the  same.  In 
practice,  a  LL  tolerance  smaller  than  the  BL  tolerance  is  chosen.  However, 
this  choice,  which  presently  worked  in  the  application  of  the  generalized 
iterative  method,  lacks  a  rigorous  rationale.  An  investigation  of  how  each 
of  the  tolerances  should  be  set  would  give  further  insight  in  this  area. 

There  is  not  a  one-to-one  correspondence  between  the  number  of 
iterations  in  the  LL  and  those  of  the  BL.  In  fact,  there  may  be  anywhere 
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from  one  to  several  hundred  Iterations  in  the  LL  for  each  iteration  in  the 
BL.  The  flowchart  in  Appendix  A  outlines  the  logic  of  this  LL-BL  scheme 
for  solving  transcendental  equations.  After  all  desired  roots  of  Eq.  (12a) 
are  found,  they  are  stored  in  an  array. 

DISPERSION  CURVES  FOR  ANTISYMMETRIC  AND  SYMMETRIC  WAVES 

The  algorithm  for  determining  roots  outlined  above  will  now  he  applied 
to  the  calculation  of  the  dispersion  curves  [3]  of  the  exact  theory,  which 
can  be  suitably  rewritten  for  the  purpose  of  these  calculations  in  the 
following  dimensionless  forms. 


For  antisymmetries!  modes,  the  equation  for  the  dispersion  curves  is 


(1  ~  aY2)(l  -  y2)**  _  (1  -  %y2)2 

tanhPsxCl  -  ay2)*5!  tanh{*s(l  -  y2)*5] 


(15) 


which  holds  when  0  <  y2  <  1,  and  for  symmetric  modes  the  equation  Is 

(1  -  av2)h(\2  -  1)**  ,  (1  -  *SY2)2 _ 

tanfexCy2  -  l)*5}  tanh[%x(l  -  ay2)*5} 


(16) 


which  holds  when  1  <  y2  <  JL  »  an^ 

a 

(1  -  *4Y2)2 _ 

tanh[Jsx(l  ~  ay2)*5 


(1  -  ay2)^!  ~  Y2)** 
tanh [Jjx(l  -  y2)*5? 


(17) 


which  holds  when  y2  <  1*  In  Eqs. (15)  through  (17),  the  normalized  wave  speed  y 
is  a  function  of  x  *  nfd/c,  where  c,  f,  and  d  are,  respectively,  the  wave 
speed  in  the  plate,  the  frequency,  and  one-half  the  thickness  of  the  plate. 
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In  Eqs.  (15)  through  (17),  a  Is  expressed  in  terms  of  Poisson's  ratio  v  of 
the  plate  material  by  the  equation 

q  -  JL-?y-  08) 

2(1  -  v)  • 

It  should  be  noted  that  In  the  dispersion  equations  for  the  lowest  order 
symmetric  Lamb  mode  [(Eqs, (16)  and  (17)],  the  value  y2  never  exceeds  1/a 
since  0  £  v  _<  0.5.  However,  in  higher  modes,  y2  does  exceed  2.0,  thus  making 
1/a  an  important  consideration. 

A  computer  program  was  designed  so  that  It  could  easily  be  implemented 
for  various  equations  with  a  subprogram  containing  the  functions  and  their 
inverses  independent  of  the  LL-BL  iterative  scheme.  In  the  subprogram,  para¬ 
meters  are  represented  by  p-numbers  (parameter  numbers).  These  p-numbers 
[4]  are  unique  to  the  computer  (HP  9825)  language  that  was  used.  This 
enables  one  to  call  the  subprogram  repeatedly  using  different  values  for  the 
parameters  each  time.  The  dispersion  curves  for  antisymmetrical  waves  are 
shown  in  Fig.  3  for  values  of  Poisson's  ratio  equal  to  0.5  (rubber),  0.303 
(steel),  and  0.05  (beryllium). 

In  order  to  calculate  these  dispersion  curves,  Eq.  (15)  must  he  broken 
down  into  four  functions  shown  in  Eq.  (12a).  The  function  assignment  chosen 
for  Eq.  (15)  is  A  *  (1  -  ay2)(l  ~  Y2)*5,  B  -  tanh[*jx(l  -  ay2)*5],  C  *  (1  -  *5Y2)2, 
and  D  ■  tanh[lj(l  -  y2)^].  The  functions  do  not  have  to  be  labeled  this  way; 
but,  whatever  choice  is  made,  A/B  must  equal  C/D.  Appendix  (D)  illustrates 
the  subprogram  for  calculating  the  dispersion  curves  associated  with  the 
first  order  antisymmetric  Lamb  mode.  The  expressions  AA,  BB,  CC,  and  DD  are 
the  inverses  of  the  corresponding  functions  A,  B,  C,  and  D.  Since  both  tanh 
and  arctanh  are  needed  for  this  set  of  functions  and  since  the  computer 
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(Hewlett-Packard  9825)  does  not  have  written  algorithms  for  these  functions, 
appropriate  subroutines  were  written  and  included  in  this  subprogram. 


*  •  0.5  (a)  r  •  0.503  *  *0.05 


Fig.  3  -  Antisymmetric  flexural  waves  of  rubber  (v  ■  0.5),  steel 
(v  ■  0.303),  and  beryllium  (v  »  0.05). 

Curve  a  is  the  flexural  wave  of  thick-plate  theory. 

Curve  b  is  the  flexural  wave  of  classical  plate  theory  for 
the  same  material. 

The  iterative  program  (see  Appendix  C)  using  the  set  of  functions  A, 

B,  C,  and  D  given  above,  performs  adequately  for  values  of  x  >  0.4.  For 
values  of  x  <  0.4,  however,  excessive  round-off  errors  were  introduced.  To 
circumvent  these  errors,  the  classical-plate-theory  equation  was  used  for 
values  of  x  <  0.4.  The  equation  for  this  dispersion  curve  is 

1  •  <19> 


This  expression  gives  a  dispersion  curve  that  is  a  straight  line,  as  shown 
by  curve  (b)  of  Fig.  3.  Curve  (a)  is  the  dispersion  curve  given  by  the 


exact  theory.  Thus,  as  Fig.  3  illustrates,  Eq.  (19)  is  a  good  approximation 
of  the  thick-plate-theory  for  x  ^  0.4.  When  calculating  a  dispersion  curve 
for  antisymmetric  waves,  the  old  root  was  used  as  the  new  seed. 

Each  dispersion  curve  for  symmetrical  waves  (see  Fig.  4)  requires  two 
sets  of  functions.  One  set  of  functions  is  entered  when  1  <  y2  <  1/a  and 
another  set  of  functions  when  y2  <  1.  Equations  (16)  and  (17)  are  the 
characteristic  equations  for  these  dispersion  curves.  To  calculate  the 
curves  shown  in  Fig.  4,  it  was  necessary  to  overcome  several  difficulties. 


DIMENSIONLESS  WANE  NUMBER  (X) 

Fig.  4.  Symmetrical  compression  waves  of  rubber  (v  ■  0.5), 
steel  (v  ■  0.303),  and  beryllium  (v  *  0.05). 

The  line  v  •  1  is  a  critical  value. 


This  was  done  by  using  a  seed  that  is  a  better  approximation  of  the  true 
root  z;  rather  than  using  a  root  from  the  previous  x  as  n  seed  for  current 
x  as  was  done  in  antisymmetric  case  [see  the  example  illustrated  using  F.q. 
(9)].  By  linearly  extrapolating  from  the  old  roots,  one  could  calculate  a 


a  root  y  for  any  x  using  the  Iterative  rootfinder.  To  determine  a  seed  Yq 
for  current  calculations  of  Yn*  one  must  compute  Ay*  The  expression  Ay  is 
given  as 


where  Y  .  is  the  root  associated  with  the  previous  valued  x  and  y  _  is 
n~i  n“<£ 

that  root  calculated  before  Yn_j*  To  compute  Yq  for  the  symmetric  mode, 
the  following  equation  is  used, 


The  algorithm  for  performing  the  extrapolation  can  be  found  in  Appendix  C 
in  the  subroutine  'DO'.  It  has  been  shown  that  as  x  increases,  'DO' 
appropriately  adjusts  f(x)  so  that  the  desired  root  is  calculated.  The 
several  difficulties  that  were  overcome  by  the  extrapolation  will  be  ex¬ 


plained  in  detail. 


The  first  difficulty  arose  when  it  was  necessary  to  change  from  one 
functional  form  to  another  at  Y  ■  1  (see  Fig.  4).  The  dispersion  curves  for 
the  symmetrical  waves  are  calculated  by  using  Eq.  (16)  when  y  <  1  and  using 
Eq.  (17)  when  1  <  y  <  a  •  However,  if  a  value  of  y  that  is  greater  than 
1  is  used  for  the  seed  of  the  iterative  scheme  at  an  x  value  where  the  true 
Y  should  be  less  than  1,  the  computer  program  will  always  work  with  the 
particular  functions  that  appear  as  A,  B,  C,  and  D  in  Eq.  (16).  These 
functions  will  always  produce  roots  Y  greater  than  1.  The  iterative  program, 
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therefore,  will  never  get  to  use  the  set  of  functions  that  appear  in  Eq. 

(17),  which  are  needed  to  obtain  a  value  of  y  less  than  1.  So  the  itera¬ 
tion  cannot  be  started  with  a  seed  y  >  1  if  it  is  to  produce  a  root  y  <  1. 

The  extrapolation  remedied  this  problem.  It  reduced  the  previously  obtained 
value  of  Y  enough  so  that  the  seed  needed  to  compute  values  of  Y  <  1,  such 
that  Y  would  itself  be  less  than  1  at  the  appropriate  x,  and  hence  allowed 
the  new  set  of  functions  to  be  entered. 

Another  difficulty  was  revealed  when  the  absolute  value  of  the  slopes  of 
A(x)  C(x) 

BOO  arU*  dTx)  were  calcu^ated  using  the  root  obtained  for  the  previous  value 
of  x  as  the  seed  for  the  calculation  with  an  incremented  x.  If  the  deriva¬ 
tives  of  A/B  and  C/D  are  nearly  equal  in  magnitude,  the  value  chosen  for  the 
seed  could  make  a  significant  difference  as  to  which  side  of  Eq.  (12a)  is 
fixed.  If  the  wrong  set  of  functions  [(i.e.,  the  wrong  side  of  Eq.  (12a)]  is 
fixed  in  the  BL,  then  the  LL  calculations  cannot  converge  to  the  true  root, 
and  the  LL  iterations  will  continue  indefinitely.  For  example,  if  Eq.  (13a) 
were  incorrectly  chosen  by  the  program  instead  of  Eq.  (13b)  in  the  BL,  the  LL 
would  be  a  root  of  Eq.  (14a)  and  (14b)  without  converging  to  a  definite  value 
that  would  return  the  calculation  to  the  BL.  However,  if  the  extrapolation 
procedure  is  used  to  determine  che  seed,  the  BL  slope  comparison  is  more 
likely  to  identify  the  correct  side  of  Eq.  (12a)  to  fix,  since  the  seed 
obtained  by  extrapolation  is  closer  to  the  true  root  than  is  the  value  of  y 
obtained  on  the  previous  calculation. 

And  the  third  problem  was  exhibited  at  values  of  x  where  the  dispersion 
curve  expressed  by  Eq.  (16)  and  the  curve  given  by 

Y?'  =  +  (ir/x) 2  (20) 
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pass  close  to  one  another.  Equation  (20)  gives  the  locus  of  values  for 
which  the  argument  of  the  tangent  in  the  equation 

B(x)  -  tanfex(y2  -  l)*5]  (21) 

is  equal  to  tt/2;  that  is,  it  gives  values  of  y  for  which  J$x(y^  -  1)  ■  tt/2. 

The  tangent  is  not  continuous  at  the  point  x  *  tt/2;  one  has 

litn  tan(t)  =  -oo  ,  and  lim  tan(t)  =  +  *> 


t  «  Ssx(y2  -  l)h  . 


It  is  important  not  to  allow  the  iterative  calculation  procedure  to  en¬ 
counter  this  point  of  discontinuity.  For  example,  suppose  x  has  been  pro¬ 
gressively  incremented  and  many  roots  have  been  calculated  from  Eq.  (16) 
with  Y2  <  1  +(^)2  . 

If  the  next  increment  of  x  is  such  that 

Y2  >  1  +  <”>'  > 

for  the  seed  used,  then  a  root  Yq  will  be  produced  that  involves  a  different 
branch  of  the  inverse  tangent  function  than  that  which  was  used  in  the 
calculation  of  previous  points  on  the  dispersion  curve.  If  this  occurs, 

the  computing  algorithm  will  produce  an  undeslred  root.  The  extrapolation 

2 

procedure,  however,  keeps  y2  <  1  +  (^)  by  subtracting  an  appropriate 
amount  from  the  Y  associated  with  the  previous  value  of  x  and  one  always 
uses  the  correct  branch  of  the  Inverse  tangent  function. 
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It  has  been  shown  all  the  difficulties  encountered  in  determining  the 
dispersion  curves  for  symmetric  waves  were  overcome  by  linearly  extra¬ 
polating.  The  iterative  program  works  for  a  wide  range  of  values  of 
Poisson's  ratio;  however,  for  small  values  of  this  quantity  (v  _<  0.05) 
additional  program  modifications  were  necessary.  These  modifications  in¬ 
clude  choosing  a  seed  for  each  x  that  is  less  than  1  +  (n/x)2.  This  was 
accomplished  by  taking  the  extrapolated  seed  e  and  checking  if 
e  <  1  +  (n/x)“.  If  so,  e  is  returned  as  the  seed  for  the  new  iteration. 
However,  if  e  >  1  +  (n/x)2,  the  quantity  1  +  (n/x)2  will  become  the  new 
seed  instead  of  e.  This  subroutine  is  called  "TN"  and  may  be  found  in 
Appendix  E.  Also,  if  the  similar-slope  problem  previously  encountered  is 
not  resolved  by  the  extrapolat ion  described  above,  further  control  is  needed 
to  determine  which  set  of  functions  in  the  BL  should  be  fixed.  This  was 
accomplished  by  fixing  one  set  of  functions  in  the  BL  and  iterate  in  the  LL 
until  a  root  was  outputted.  However,  if  a  divergence  exists  in  successive 
approximations,  the  wrong  set  of  functions  was  fixed  in  the  BL.  Therefore, 
the  other  set  of  functions  would  be  fixed.  A  divergence  is  assumed  when 
the  number  of  LL  iterations  exceeds  500.  The  number  500  is  arbitrary,  and 
possibly  a  larger  number  would  be  needed  for  some  value  of  Poisson's  ratio. 
These  changes  are  also  found  in  the  program  in  Appendix  E. 

This  revised  program  worked  for  beryllium  (v  =  0.05)  and  is  believed 
to  work  for  v  <  0.05.  This  alternative  iterative  program  can  be  found  in 
Appendix  E  as  a  supplement  to  the  original  program  in  Appendix  C.  However, 
the  modification  of  the  program  produces  the  roots  much  more  slowly.  For 
this  reason,  it  is  recommended  that  one  use  the  program  in  Appendix  C 
unless  they  are  using  a  value  of  Poisson's  ratio  less  than  0.10. 
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ADDITONAL  CALCULATIONS 

Kappas  (<) 

Figure  5  illustrates  dispersion  curves  for  both  the  flexural  and  compres- 
sional  waves  of  a  rubber  material  (v  *  0.5).  Both  curves  approach  a  real 
value  =>  cR/cg  where  c^  and  cg  are,  respectively,  the  Rayleigh  wave  speed 
and  the  shear  wave  speed. 

2 

It  was  of  interest  to  graphically  illustrate  the  behavior  of  <  as  a 
function  of  Poisson's  ratio  when  k  is  computed  by  two  different  methods. 


DIMENSIONLESS  WA/E  NUMBER  (X) 


Fig.  5.  Symmetrical  [1]  and  antisymmetrical  [2]  waves  of  rubber. 

Kappa  is  the  point  where  the  two  curves  approach  the 
normalized  Rayleigh  wave  velocity 


Note  chat  as  ex  changes  from  0  to  0.5,  Poisson's  ratio  v  changes  from  0.5  to 
0.  Exact  elasticity  theory,  ic  gives  as  a  real  root  of  the  equation 

<6  -  8k4  +  8(3  -  2ot)ic?  -  16(1  -  a)  =  0 

where  a  is  defined  in  Eq.  (18).  The  real  root  of  this  equation  is  plotted 
as  the  solid  curve  in  Fig.  6.  In  the  manuafacturer-supplied  software  package 
for  the  HP  9825  computer,  there  was  a  "Polynomial  Rootfinder"  program  that 


Fig.  6.  Viktorov's  approximation  of  <2  compared  with 

the  exact  theory  value  for  k2  as  a  function  of  a. 


could  determine  the  roots  of  such  polynomial  equations.  However,  in  that 
program,  the  coefficients  of  the  polynomial  must  be  manually  entered  and 
the  output  will  consist  of  all  the  roots  (both  real  and  complex).  To 
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calculate  k  as  a  function  of  the  variable  a,  It  was  necessary  to  have  the 
computer  determine  the  coefficients  and  choose  only  the  real  roots. 

The  approximation  of  k  proposed  by  Viktorov  [5]  is  given  by 


By  calculating  K  as  a  function  of  a  using  the  results  of  Eq.  (18),  one 
obtains  the  dashed  curve  in  Fig.  6.  Note  the  agreement  between  the  two 
curves  in  Fig.  6. 

Inflection  Points 

The  quantity  k  has  an  important  role  in  approximate  plate  theory. 

Dr.  P.  Dubbelday  is  investigating  improved  expressions  for  k.  In  this 
investigation,  it  was  necessary  to  examine  features  of  the  strain  distri¬ 
bution  for  antisymmetric  waves.  The  formula  for  the  strain  distribution 
produced  by  antisymmetric  waves  is 

ezx  =  cosh[qd(z/d) lcosh(sd)  -  cosh[sd(z/d)]cosh(qd)  (22) 

where  z  is  the  distance  from  center  of  the  plate  and  d  is  half  the  plate 
thickness.  The  expressions  qd  and  sd  are  given  by  the  equations  qd  = 
x(l  -  ay2)'5  and  sd  =  x(l  -  y2)'5  whereas  previously  x  =  ( Trd ) / A  with  X  the  wave¬ 
length  of  the  antisymmetric  wave.  In  particular,  it  was  believed  that  the 
inflection  point  of  Eq.  (22)  should  be  an  identifiable  point  of  the  curve  that 
merges  with  the  center  of  the  plate  at  the  low-frequency  limit.  When  the  de¬ 
rivative  c’zx  is  equal  to  zero,  we  obtain  a  maximum  or  minimum;  and  when  the 
second  derivative  e"  is  equal  to  zero,  we  get  an  inflection  point.  The  in- 

flection  points  (kz^  are  given  in  the  form  (qd) 2cosh[kzi (l  -  ay2)J5]cosh(sd)  = 

2  2  ^ 

(sd)  coshtkz^d  -  aY  )  ]cosh(qd).  Figure  7  illustrates  the  location  of  the 
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Fig.  7.  Point  of  inflection  for  rubber  (0.5),  steel 
(0.302774),  and  beryllium  (0.05). 


inflection  point  in  the  plate  with  respect  to  the  plate’s  center  as  a  func¬ 
tion  of  x  for  three  materials,  steel  (v  =  0.302774),  rubber  (v  =  0.5),  and 
beryllium  (v  =  0.05).  It  was  originally  believed  that  these  inflection 
points  were  a  characteristic  feature  that  could  be  used  in  the  theory  being 
developed  for  calculating  «.  However,  the  graph  reveals  the  fact  that  points 
of  inflection  do  not  exist  for  values  of  x  less  than  approximately  3.2. 

Thus,  the  hypothesized  method  of  finding  k  could  not  be  used. 


Maximum  Points 

Finding  the  point  of  maximum  strain  kz  parallels  the  findings  of  in- 

m 

flection  points.  The  maximum  points  kz  are  given  as  the  root  of 


(qd)sinhfkz  (1  -  ay2  )h]  =  (:;d)sinh[kz  (1  -  Y2)!2]. 
in  m 
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Figure  8  illustrates  the  location  of  points  of  maximum  for  rubber  and 
beryllium. 


10.0 


THICKNESS  PARAMETER  (X) 

Fig.  8.  Point  of  maximum  for  rubber  (0.5)  and  beryllium  (0.05). 

Note  that  kz  is  zero  for  values  of  x  less  than  approx.  3.2. 
ra 

The  point  of  maximum  and  the  point  of  inflection  intercept  the  x  axis 
at  approximately  the  same  point  for  each  material;  however,  unlike  the  in¬ 
flection  curve,  the  values  of  the  maximum  curve  exists  for  all  x.  They  are 
zero  for  x  less  than  3.2. 

Strain  Distribution 

After  performing  the  calculations  to  find  points  of  the  inflection  and 

maximum  in  the  strain  component  that  is  given  by  Eq.  (22);  the  normalized 

strain  distribution  e  /z  ,  .  for  distinct  values  of  x  was  calculated. 

zx  zx(m) 
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The  formula  for  this  normalized  strain  distribution  is  E  /e  ,  »,  where 

zx  zx(m) 


€  ,  •*  =  cosh(qz  )cosh(sd)  -  cosh(sz  )cosh(qd) 

zxim;  ra  ra 


and  e  is  given  in  Eq.  (21).  The  value  qz  and  sz  are  given  by  the 
zx  mm 

following  equations,  qz  =  kz  (1  -  ay2)'1  and  sz  =  kz  (1  -  y2)^.  Figure  9 

mm  mm'6 


Fig.  9.  Strain  distribution  of  the  exact  theory  with  a  =  0.302774. 

illustrates  the  strain  distribution  for  various  values  of  x.  Note  that  this 
quantity  has  been  plotted  so  that  z/d  is  the  ordinate  and  the  strain  distri¬ 
bution  is  the  abscissa.  It  was  done  in  this  manner  to  best  reflect  the 
physical  situation  by  a  propagating  antisymmetric  Lamb  wave.  Note  when  x 
is  less  than  about  3.2,  the  transverse  shear  strain  is  1  in  a  region  near 
the  center  of  the  plate.  As  x  increases  above  3.2,  the  shear  strain  in  this 
region  decreases  for  values  and  the  position  of  maximum  strain  moves  towards 
the  surface  of  the  plate. 
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The  strain  distribution  is  plotted  in  Fig.  10  for  three  materials  when 
x  has  a  value  of  10.  Figure  10  shows  that  the  normalized  strain  is  quite 
different  at  the  center  of  the  plate  and  yet  much  the  same  near  the  surface. 
Similar  results  were  found  when  kd  assumed  values  other  than  10. 


TRANSVERSE  SHEAR  STRAIN 

Fig-  10.  Comparison  of  the  strain  distribution  of  three  materials 
with  x  *  10. 

SUMMARY 

The  roots  of  transcendental  equations  that  arise  in  elastic  plate  theory 
were  calculated  using  a  generalization  of  a  well  known  iterative  method.  The 
basic  iterative  method  is  successively  applied  to  combinations  of  functions 
appearing  in  the  transcendental  equations  until  the  root  with  the  desired 
tolerance  is  determined.  The  generalized  rootfinder  described  is  a  straight¬ 
forward  technique  for  solving  transcendental  or  algebraic  equations.  This 


generalized  scheme  was  applied  to  the  exact  elasticity  theory  dispersion 
equation,  the  polynomial  equations  for  calculating  and  the  equations 
that  describe  points  of  strain  inflection  and  strain  maximum.  It  was  found 
that  inflection  points  in  the  strain  distribution  do  not  exist  at  all 
frequencies.  Since  the  time  this  paper  was  drafted,  several  other  research 
projects  at  USRD  have  implemented  this  generalized  method  in  their  numerical 
analysis. 

For  this  reason,  it  appears  that  other  investigators,  particularly  those 
working  in  elasticity  theory  and  in  the  area  of  numerical  analysis  should  be 
able  to  make  use  of  the  iterative  rootfinding  method  with  minimal  difficulty. 
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APPENDIX  A 


TOLERANCE  AND  ERROR  SIZE 


It  is  noted  that  the  difference  between  two  successive  approximate 
roots  of  an  equation  is  not  necessarily  a  measure  of  the  error  between  the 
approximate  and  the  true  root  of  the  equation.  The  tolerance  will  be  defined 
as  the  difference  between  successive  approximations  and  the  error  size  de¬ 
fined  as  the  difference  between  the  approximate  and  true  root.  The  itera¬ 
tive  method  described  in  the  text  is  used  to  solve  an  equation  in  the  form 

f(x)  =  g(x).  (Al) 


It  will  be  shown  that  the  tolerance  and  error  size  can  vary  significantly 
when  the  absolute  values  of  the  slopes  of  f(x)  and  g(x)  in  Eq.  (Al)  are 
similar.  Suppose  two  lines  to  represent  curves  are  chosen  with  slopes  m^ 
and  m2  such  that 


and 


y  =  mxx  (A2) 

y  =  m2x.  (A3) 


Take  m2  less  than  m^  and  suppose  both  slopes  have  the  same  sign.  Let  z  be 
the  actual  root  of 

m^x  =  m2x.  (A4) 


Take  x^  to  be  initial  approximation  of  z.  If  x  is  replaced  by  x^  in  the 
right-hand  side  of  Eq.  (A4) ,  then 


in  jX  —  m  ^  q  j 


(A5) 
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which  can  be  solved  by  calculating  x,  in  the 


x^  =*  n^x^/m^  (A6) 

At  this  stage  the  tolerance  can  be  written  as  | x^  -  x^ |  and  the  error  size 
as  | x^  -  z|.  It  is  desirable  to  Investigate  the  situation  in  which  toler¬ 
ance  is  less  than  a  small  number  e,  and  the  error  size  is  greater  than  e. 
Symbolically,  the  following  two  relations  are  to  be  satisfied: 


1*1  ■  xol  =  lra2  x0/ml  "  xo! 


(A7) 


and 


IX1  ~  z 


|n>2  “  0 1  >  F  * 


(A8) 


Note  that  z  is  equal  to  zero  in  Eq.  (A8),  since  the  curves  expressed  in  F.qs. 
(A2)  and  (A3)  intersect  at  the  origin  [x  =  0  satisfies  Eq.  (A4)]. 

Eq.  (A7)  may  be  re-written  as 

K  ‘  V  ■  i*o  i  ^  - 11  ' E 


or 


Ixf  ~  xQ  |  =  |*0|  (!-—)<£  • 


Therefore 


m2/m1  >  1  -  e/ | Xq | . 


(A9) 


By  substituting  1  -  e/)xg|  for  n^/m^  in  Eq.  (A8),  one  obtains 


lxi  '  zi  >  U  5  >xo'  >  f  ’ 


which  may  be  simplified  to 


>  £  * 
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Hence 


IX1  "  zl  >  lxol  >  2c. 

Therefore,  If  the  seed  Xq  is  chosen  such  that  Xq  >  2e,  then  the  calculated 
root  x^  will  differ  from  Xq  by  less  than  the  tolerance  £  without  x  being 
within  a  distance  e  of  the  true  root 
greater  than  2g,  the  tolerance  |x^  - 

Take,  for  example,  the  lines 

y 

and 

y  ■ 

It  is  clear  that  the  slopes  of  these 
point  of  intersection  is  (0,0).  The 

Suppose  a  guess  of  0.5  was  made 
guess  be  called  Xq.  If  Xq  is  put  in 
If  this  value  is  put  in  the  place  of 
mation  x^  is  equal  to  0.4995. 

It  can  be  seen  that  the  difference  |Xq  ~  x^ |  ■  0.0005.  Certainly  if  e 
were  set  at  1  *  10  ,  the  difference  |xq  -  x^ |  <  e. 

However,  the  difference  |z  -  Xq|  ■  0.4995  which  is  not  less  than  the 
predetermined  tolerance.  Hence,  a  rather  simple  example  shows  why  tolerance 
is  not  always  a  good  measure  of  error. 

The  author  is  aware  that  an  argument  has  been  presented  only  for  linear 
equations  at  the  origin.  Linear  equations  were  chosen  because  curves  appear 


z.  Thus,  wherever  the  error  size  is 
Xq  |  is  not  a  good  measure  of  the  error. 


(All) 


999 

1000 


(A12) 


two  lines  are  very  similar  and  that  the 

999 


actual  root  z  of  x 


1000 


x  is  0. 


for  the  value  of  this  root.  Let  this 
Eq.  (A12),  y  takes  the  value  0.4995. 
y  in  Eq.  (All),  the  first  approxi- 


much  like  straight  lines  near  the  point  of  intersection.  It  is  chosen  that 
the  equations  have  y-intercepts  of  zero  in  the  examples  for  simplicity.  A 
parallel  argument  holds  for  curves  that  intersect  at  a  point  other  than  the 
origin. 


32 


APPENDIX  B 


GENERAL  FLOW  DIAGRAM 
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INPUT  SEQUENCE 


LITTLE  LOOP  TOLERANCE  DECISION 
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BIG  LOOP 

TOLERANCE  DECISION 


COMPLETION 

AND 

CONTINUATION 
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appendix  c 


ITERATIVE  HP  PROGRAM 

0 :  fit  9 
1 :  Q- K j  0  +  C 
2:  .  302774*  V  5 

P r  t  " V "  J  v  COMMENTS 

3:  •' ALPHA’  (V)  *r2 

4!  1  £  -  4  *J  Line  2:  Poisson's  ratio  (V) 

5:  le-6*Z 

£ !  ent  "nornO.  liz  Line  3:  Alpha  m  terms  of  Poisson's 

ed  wave  speed"  >  ratio.  See  line  47. 


rl 

7:  1 e - 5  *  R 

Line  4: 

Big  Loop  Tolerance  (J) 

8:  1  e - 5  *  M 

Little  Loop  Tolerance  (Z) 

9 :  r  1  *  r  1  *  r  1  j  d  i « 

A  Cl 00 j  2] 

Line  6: 

Original  Seed  (rl) 

10:  ent  "ent  1 

Lines 

i f  asvnroet r i c " » 

A 

7  &  8: 

Big  Loop  and  Little  Loop  deltas 

11:  for  X = 1  to 

Lines 

100  5  if  A  =  1 5  =»  =•  b 

11-14: 

Choosing  the  correct  set  of 

1 5 

functions  to  be  entered. 

12:  if  r 1 < 1 5  2^ A 5 
?sb  15 

Line  17 : 

Initializing  and  incrementing 

13:  if  r 1 > 1  a  n  d 

parameter  (r3) 

r  1  <  1  /  r  2  5  3  *  A  5 
•?sb  15 

Lines 

14:  if  r 1 > 1 / r  2  5 

18  &  19: 

Determining  the  slopes  of 

4*R  j  ?sb  15  A/B  and  C/D 

15s  if  A  =  C  j  • j  n  p  2 

16:  R*c;idf  A » 

5 1  j  17 

17:  "He  re":. 2  + 

.2X*r3Jprt  "X"» 
r3 

18:  a  b  s ( ’ A ’ ( r  1  + 

Rj  r  2  j  r  3 )  /  ’  B  ’  (  rl 
+  R> r 2r r  3 ) - ’ A ’  (  r 
1 -R  j r 2  j r3 ) / 5  B J  ( 
r  1  -  R  j  r  2  j  r  3 )  )  *  U 
19:  o.bs  (’O’  (  rl  + 

R,  r2>  r 3 )  / ’  D ’  (  rl 
+  R  >  r2  >  r  3 )  -  ■'  C»  (  r 
1 -R  j r 2  » r 3 ) /  ’  D ’  ( 
r  1  -  R  j  r  2  ?  r  3 )  )  *  T 
20:  rl*S5r3*Y 
21 :  if  U - T <  =  0 5 
■  j  n  p  2 

22:  if  T - U < 0 > 
jfiP  8 
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4 


23:  »RUS»  r2» 

42:  if  abs  (EX*  J 

Y)  +r8i  »  BMS»  r 2  j 

5  9 1  o  "print" 

Y )  *  r  9  5  r  3  /  r  9  *  F 

43:  I*r  lXsb  13 

24:  absU*DMS  + 

44:  "print" : r3*A 

M ? r 2  j  Y ) - ’ 0 ’  (S- 

ex ?  n ;  thacxi 

M» r 2  j  Y ) ) *F> +  0 

23  5  P  r  t  "  R  "  »  f  1 5 

25:  abs  ( !  CMS  +  Mi 

I-’ DQMI  Hr  15 

r 2 j Y) -»C»  (S-Mi 

0  +  U5  0-»Ki  rl  +  r  135 

r  2  j  Y )  )  *P 

n  e  x  t  X 

26:  if  O-P<=05 

45:  ent.  "  ldf 

•  j  n  p  2 

#for  Data"  <  B 

27:  j  m  p.  2 

46:  trk  0  5  ref  B> 

23:  -1 1  ->L  5  J  O’  (Si 

fi  C  •*  3  5  2nd 

r 2s Y) *F*Q5 > CC’  ( 

47:  "ALPHA" s ret 

3  ?  r2iYiQHI5 

(1-2V) /2(1-V) 

9 t o  "LL" 

43:  "DO":  Hr 5 

29:  -104L 5 > C*  (St 

49:  if  X= 1  or 

r 2  j  Y ) /F*Q5 » DO’  ( 

A  =  1 5  r  5  *  r  6  5  r  e  t  0 

Sj r 2  j  Y  i  Q  H 1 5 

50:  r6-r5+r75 

9 to  "LL" 

r  5  *  r  6  5  r  *  t  r  7 

30:  »CMSir2iY)/ 

1  DM3)  r2i  Y)  +  F 

31:  abstt’B’  «S+ 

Mi r 2  j  Y ) - ’ B ’  (3- 

Mi  r2>  Y)  )  *  F )  *  0 
32:  abs(»flMS  +  M» 
r  2  j  Y )  -  ’  R  ’  ( 3  -  M  j 
r  2  ?  Y )  HP 
33:  if  0 - P >  =  O 5 
j  ivi  p  2 


Lines 


COMMENTS 


23  &  30:  Fixing  equation  with  smaller  slope 
Lines  24,  25, 

31,  &  32:  Slope  Decision  in  Little  Loop 


34:  -5*L5JBMSi 
r2»  Y)  *F*Q!  ’RAM 
3  j  r  2  j  Y  >  Q  H 1 5 
9to  "LLM 
35:  -4*U  *  fl»  (S» 
r  2  >  Y )  /  F  *  Q  5 ’ BB ’  ( 
3  j  r  2 »  Y  >  Q )  ■*  1 5 

9 1  o  "LL" 

36:  "LL": 

37:  S-HG 

33:  if  abs (G) <=2 
5K+HK!*to  "BL" 
39:  H  3  5 . j  i  vi  p  L 
40:  "  B  L " :  if  K= 1 5 
r  1  r  1 4 5  Hr  1 5 

92  b  13 
41:  rl-HE 


Lines 

26  &  33:  LL  Slope  comparison 

Lines  28,  29, 

34,  &  35:  Solving  for  variable  with  smaller 
slope 

Line  38:  Little  Loop  Tolerance  Decision 

Line  42:  Big  Loop  Tolerance  Decision 

Line  44:  Prints  roots  (rl)  and  incremental 

parameter  (r3).  The  subroutine 
’DO’  is  a  linear  interpolation  for 
new  seed.  See  Lines  48-50. 

Lines 

45  &  46:  Loads  and  stores  array  on  tape 


40 


p9 

9 :  "  B  B " :  ’  fi  r  t  a  n  h  * 
(  p  4  )  / .  5  p  3  *  p  5 
10:  1-p5I'2*p6 

11:  Pb / p 2  * p  7  5 
r  ■=•  t,  p  7 

12:  "CC" : if  p1>2 
j  •  j  ivi  p  2 

13:  2*(l-rabs(p4 
J  )  *  p  5  5  r  e  X  p5 
14:  2* ( 1+rabs (p4 
1  )  4  p  5  5  r  e  t-  p  5 
15:  "DD" : ’ firtcmh 
’  (  p  4  J  /  ( .  5  *  p  3 )  -*  p 
5 


1  6 : 

p 

5 1 2  *  p  6 

17: 

r. 

i 

-  p  6  ■*  p  7 

• 

J 

r 

e  t. 

P'  1 
18: 

•• 

T  a  n  h " : 

( 

1 

- 

*E‘  X 

p 

( 

-  2  p  1 )  ) 

i 

1  + 

€'  X 

p 

t 

-2p1  J  ) 

-♦ 

P 

9  5 

re 

t 

p9 

19: 

•I 

R  r  t.  a  n  h 

M 

■ 

■ 

f  1  + 

p  1 

;i 

( 1-Pl ) 

p 

8 

2  0 : 

■ 

5*ln(P 

8 

1 

+  p9 

j  r 

% 

p9 

*13 

3 

5 

3 
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APPENDIX  E 


ITERATIVE  HP  PROGRAM  FOR  v  <  10 


0  ■  f  1 1-  9 
1 :  0  ■>  K  j  0  -*■  C 
2 •  .05* V ! p r t 

"  V  "  j  V 

3:  J ALPHA’ (V) *r2 
4  :  1  e  -  6  *  J 

5:  le-6*Z 
6  i  €'  n  t  “  n  o  r  n  o.  1  i  z 
e  d  i.'j  o.  v  €•  2-  p  6 1  d "  j 


r  1 

7:  U-5  +  R 
3:  U-5*M 
9  •  r  1  r  1  ->  r  1 5  d  i  f*i 
fl  C 1  0  0 »  2  3 
10:  e n t  " e  n  t  1 
i  f  o.  2  y  m  n  *■ t-  r  i  c "  j 
fl 

11:  f r  X  =  1  t  o 
1  00  5  if  fl  =  1 5  9 2- b 
15 

12:  it  r  1  <  1 5  2 * ft 5 
=>sb  15 

13:  i f  r 1 > 1  a n d 
r  1  <  1  /  r  2  J  3  *  fl  ? 

•ssb  15 

14:  if  r 1 > 1 / r  2  5 
4  •>  R  j  9  2-  b  15 
15:  if  ft  =  C 5  -i ft p  2 
16:  fl*C51df  fl? 


.  3  a  •*  r  i  j  p  r  l  a  j 
r 3 5  ’  TN’  (  rl  j  r3)  * 
r  1 

13:  abs ( * fl* ( r  1  + 

R  j  r  2  >  r  3 )  /  ’  B  ■  (  r  1 
+  R  ?  r  2  ?  r  3 )  -  ’  fi  ’  (  r 
1-R  j  r2  »  r 3 )  5  6  •  ( 

r  1  -  R  ?  r  2  j  r  3 )  )  *  U 
19:  abs ( ’ C ’  (  r  1  + 

R  j  r  2  j  r  3 )  /  ’  D  ’  •  r  1 
+  R  ?  r  2  j  r  3 )  -  ’  C  ’  (  r 
1  -  R  ? r 2  ? r  3 ) - ’ D ’  ( 
r  1  -  R »  r  2  ?  r  3 )  j  *  T 
20:  r 1*3? r  3  * Y 
21 :  if  U - T <  =  0 5 
l*N5-17+r315 
j  ft  p  2 


COMMENTS 


Line  17;  To  avoid  interfering  with  tan 
curve,  'TN'  was  designed  to 
adjust  resubmitted  seed. 

See  Lines  51-53. 
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.  -vyg? 


• 

w  • 

if  T-U 

<0 

a 

1 

23N 

j  -  1 0  ■*  r 

31 

5 

•  j  ivi  p 

r! 

■“«  a 
•— 1  • 

’  ft  ’  ( s » 

r2 

5 

Y  j  3 

r  3  5  ’  B  ' 

i  s 

j  r2 

Y )  - 

r  9 f  y 

r  9 

3F 

4: 

o.  b  s  (  ( ’ 

D  * 

( 8  + 

M  j  r 

2  j  Y  :i  -  » 

D 9 

(8- 

M  j  r 

2  ?  y  :i )  * 

f:i 

3  0 

5: 

abs  C 5  C 

» ( 

S  +  M 

r  2  j 

Y )  - »  C  ’ 

(8 

-M  j 

r  2  > 

Y)  )  3P 

6 : 

i  f  0  -  P 

.(  — 

0  5 

j  o  p 

•-i 

C. 

■?  • 

l  • 

o  i  vi  p  2 

ij  a 
'•*  * 

-1 1+LS 

J  D 

»  (S 

r2  j 

Y) #F3Q 

; » 

cc» 

8 1  r 

2  j  Y  ?  Q  J 

■»I 

a 

9 

3  t  0 

"LL" 

9: 

-103LS 

*  c 

5  (S 

r2  > 

Y )  F  3  Q 

; » 

DO* 

8  j  r 

2  j  Y  j  Q ) 

4  I 

■ 

J 

3 1  0 

"LL" 

0 : 

>  C  ■  ( 8  j 

r2 

j  Y ) 

’O’ 

( S  j  r  2  ? 

Y ) 

3F 

1 : 

abs  (  ( ’ 

B  ’ 

(8  + 

M  j  r 

2  j  Y )  - » 

B» 

i:s- 

M  j  r 

2  j  Y )  j  * 

F) 

3  0 

•**i  a 

w  ■ 

abs  ( *  ft 

»  ( 

S  +  M 

r  2  j 

Y )  -  »  ft  ’ 

i;  s 

-M> 

r2  ? 

Y)  )  3p 

3: 

if  0-P 

■>  = 

0 ; 

•  j  iyi  p 

o 

w 

4: 

-5^l;  j 

B» 

( 3  j 

r2  j 

Y ) *F* Q 

•  1 

J 

flft’ 

8 ) r2> Y  > Q J *If 
3to  "LL" 

35:  -4*L;jflUSj 

r 2  j  Y  J  /F3QJ  J  BE: J  ( 

S  j  r  2  >  Y  >  Q )  3  I  > 

3*0  "LL" 

36:  " L  L " : 8 - 1 +  G  5 
U+l3U»if  M>500 
o.  n  d  1  =  N  5  r  1 4  3 1 J  - 
1 0  3  r  3 1  >  0  3  U  j  3 1.  o 
43 


COMMENTS 

Line  36:  Number  of  LL  iterations  (W) 
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37:  if  U> 500 
and  2  =  N  5  r  1 4  *  1 5 
•g  t  o  43 

3*3:  if  abs  (G)  <=2 

;  0-hj;  k+i^k;  3to 

"EL" 

39:  I3S5omp  L 
4@:  “BL“ : if  K=  1 5 
rl+rl4SI*rl*S» 

J ivip  r31 
41:  r  1  —  I  ->  E 
42:  if  a b s ( E ) <  =  J 
5  3 1  o  "print" 

43:  I  ■>  r  1  ■>  S  5  J  ft  p 
r3 1  -3 

44:  "print": r3+A 
cx>  n  ;n*ACX» 

2] 5  prt  “R">ri; 
I-’DOMD+rl? 
03U!0->K;rl->rl35 
next  X 

45:  ent  "ldf 
#for  Data“>B 
46:  t  r k  0  5  ref  B> 
fl [*] j  end 
47:  "ALPHA" : ret 
(1-2V)/2(1-V) 
48:  “DO” :  I*r5 
49:  if  X  = 1  o r 
A  =  1 5  r  5  *  r  6  5  r  e  t  0 
50:  r6-r5+r7> 
r  5  *  r  6  5  r  e  t  r  7 
51:  ”TN" :  1  +  (if/ 
r  3 )  1 2  *  r  6  0 
52:  if  r68<r 1 5 
ret  r6@ 

53:  ret  r  1 
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